custom_theme <- theme_bw() + 
  theme(text = element_text(family = "Times New Roman", size = 10))
theme_set(custom_theme)

tanzania_pal <- rep(brewer.pal(9,"Pastel1"),times=4)
market_pal <- c(brewer.pal(8, "Dark2"), "#FFA500", "#6A3D9A", "#33A02C")

Tanzania Map

tregs_shp <- tregs_shp %>%
  rename(
    "Country_N"=ADM0_EN,    
    "Country_NSw" = ADM0_SW,
    "Country_C"=ADM0_PCODE,
    "Region_Nam" = ADM1_EN,
    "Region_Cod" = ADM1_PCODE,
  )

twards_shp <- twards_shp %>%
  rename(
    "Country_N"=ADM0_EN,    
    "Country_NSw" = ADM0_SW,
    "Country_C"=ADM0_PCODE,
    "Region_Nam" = ADM1_EN,
    "Region_Cod" = ADM1_PCODE,
    "District_N" = ADM2_EN,
    "District_C" = ADM2_PCODE,
    "Ward_Name"=ADM3_EN,
    "Ward_Code"= ADM3_PCODE
  )

tanzania_map <- ggplot() +
  geom_sf(data = tregs_shp, aes(fill = Region_Nam)) +
  # geom_sf(data = tanzania_shp, aes(fill = Region_Nam, col = factor(District_C))) +
  geom_sf_text(data = tregs_shp, aes(label = Region_Nam), size = 2)+
  scale_fill_manual(values=tanzania_pal, guide = "none")+
  scale_color_discrete(guide = "none")+
  labs(title = "Tanzania map of Regions")

tanzania_map
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

twards_map <- ggplot() +
  geom_sf(data = twards_shp, aes(fill = Region_Nam)) +
  scale_fill_manual(values=tanzania_pal)+  
  scale_color_discrete(guide = "none")+
  labs(title = "Tanzania map of Regions (colour) and Ward (boundary)")

twards_map

# Market Maps (1) : Manual Clean

Market A Map - Manually Cleaned.

  • Market A map survey locations data consistent with regional tanzania map for all surveys except survey in Ruvuma, labelled as Rukwa.
# convert mrkta location data to mapping format:

mrkta_mapdata <- mrkta_gps %>% 
  mutate(lat = str_split(`gps.coordinates.of.the.market`, "\\s", simplify = TRUE)[, 1],
         long = str_split(`gps.coordinates.of.the.market`, "\\s", simplify = TRUE)[, 2]) %>% 
  filter(!is.na(lat), !is.na(long))

# for now ignore rows which have na in lat or long:
mrkta_sf <- st_as_sf(mrkta_mapdata, coords = c("long","lat"),  crs = 4326)

mrktA_map <- ggplot() +
  geom_sf(data = tregs_shp, aes(fill = Region_Nam), linewidth = 0.1) +
  scale_fill_manual(values=alpha(tanzania_pal,0.7), guide = "none")+
  geom_sf(data = mrkta_sf, aes(col = district.region.name), size = 1.5)+
  geom_sf_text(data = tregs_shp, aes(label = Region_Nam), size = 2)+
  scale_colour_manual(values = market_pal)+
  coord_sf()+
  labs(title = "Market survey A: survey locations", 
       fill = "Tanzania region",
       col = "Market region")+
  theme(legend.position = "right", legend.box = "horizontal")

mrktA_map
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

Market C Map - manually cleaned.

  • Manual cleaning, market C cross-referenced with market A (not geo-matched)
mrktc_mapdata <- mrktc_gps %>% 
  mutate(lat = str_split(`gps.coordinates.of.the.market`, "\\s", simplify = TRUE)[, 1],
         long = str_split(`gps.coordinates.of.the.market`, "\\s", simplify = TRUE)[, 2]) %>% 
  arrange(gps.coordinates.of.the.market) %>% 
  filter(!is.na(lat), !is.na(long))

# for now ignore rows which have na in lat or long:
mrktc_sf <- st_as_sf(mrktc_mapdata, coords = c("long","lat"),  crs = 4326) 

mrktC_map <- ggplot() +
  geom_sf(data = tregs_shp, aes(fill = Region_Nam), linewidth = 0.1) +
  scale_fill_manual(values=alpha(tanzania_pal,0.7), guide = "none")+
  geom_sf(data = mrktc_sf, aes(col = district.region.name), size = 1)+
  geom_sf_text(data = tregs_shp, aes(label = Region_Nam), size = 2)+
  scale_colour_manual(values = market_pal)+
  coord_sf()+ 
  labs(title = "Market survey C: survey locations", 
       fill = "Tanzania region",
       col = "Market region")+
  theme(legend.position = "right", legend.box = "horizontal")

mrktC_map
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

mrktC3andA_map <- ggplot() +
   geom_sf(data = tregs_shp, aes(fill = Region_Nam), linewidth = 0.1) +
   scale_fill_manual(values=alpha(tanzania_pal,0.7), guide = "none")+
   geom_sf(data = mrkta_sf, aes(col = district.region.name), alpha = 0.5, size = 3)+
   geom_sf(data = mrktc_sf, aes(col = district.region.name),  shape = 17, size = 1.5)+
  geom_sf_text(data = tregs_shp, aes(label = Region_Nam), size = 2)+
   scale_colour_manual(values = market_pal)+
   coord_sf()+
   labs(title = "Market Survey A (circles) and C (triangles) - manually cleaned location data",
        fill = "Tanzania region",
        col = "Market region")+
  theme(legend.position = "right", legend.box = "horizontal")
 
mrktC3andA_map
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

Market C Map 2 - market C geo-matched to market A

mrktCtoAmatch <- mrktc_mapdata %>%
  mutate(mrktaID = st_nearest_feature(mrktc_sf, mrkta_sf)) %>% 
  left_join(
    mrkta_mapdata %>%
      mutate(mrktaID = 1:nrow(.)) %>%
      select(
        mrktaID,
        "district.region.code.A" = "district.region.code",
        "district.region.name.A" = "district.region.name",
        "ward.code.A" = "ward.code",
        "ward.name.A" = "ward.name" ,
        "village.A" = "village",
        "market.name.A" = "market.name.corrected"
      ),
    by = c("mrktaID"))


mrktCtoAmatch_gps <- mrktCtoAmatch %>%
  select(
      country,
      ends_with(".A"),
      lat,
      long,
      gps.coordinates.of.the.market
    )
  
mrktCtoAmatch_sf <- st_as_sf(mrktCtoAmatch_gps, coords = c("long","lat"),  crs = 4326)

mrktCAmatched_map <- ggplot() +
   geom_sf(data = tregs_shp, aes(fill = Region_Nam), linewidth = 0.1) +
   scale_fill_manual(values=alpha(tanzania_pal,0.7), guide = "none")+
   geom_sf(data = mrkta_sf, aes(col = district.region.name), alpha = 0.5, size = 3)+
   geom_sf(data = mrktCtoAmatch_sf, aes(col = district.region.name.A),  shape = 17, size = 1.5)+
  geom_sf_text(data = tregs_shp, aes(label = Region_Nam), size = 2)+
   scale_colour_manual(values = market_pal)+
   coord_sf()+
   labs(title = "Market Survey A (circles) and C (triangles) - market locations",
        fill = "Tanzania region",
        col = "Market region")+
  theme(legend.position = "right", legend.box = "horizontal")
 
mrktCAmatched_map
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

###################
# PLOTS BY REGION #
###################
# mrkta_sf <- mrkta_sf %>% mutate(mrktID = paste(ward.name,village,market.name.corrected, sep = "-"))
# mrktCtoAmatch_sf <- mrktCtoAmatch_sf %>% mutate(mrktID = paste(ward.name.A,village.A,market.name.A, sep = "-"))
# 
# districts <- mrktc_sf %>% distinct(district.region.name) %>% pull()
# Districts <- str_to_title(districts)
# 
# map_plots <- list()
# for (i in 1:length(districts)) {
#   
#   District_i <- Districts[i]
#   district_i <- districts[i]
#   
#   map_plots[[i]] <- ggplot() +
#     geom_sf(data = tregs_shp %>% filter(Region_Nam == District_i), fill = NA) +
#     geom_sf(data = mrkta_sf %>% filter(district.region.name == district_i), aes(col = mrktID), alpha = 0.5, size = 3)+
#     geom_sf(data = mrktCtoAmatch_sf %>% filter(district.region.name.A == district_i), aes(col = mrktID), shape = 17, size = 1)+
#     scale_colour_brewer(palette = "Paired")+
#     facet_wrap(~district.region.name)+
#     coord_sf()+
#     theme_bw()
# }
# 
# map_plots[[1]] <- ggplot() +
#   geom_sf(data = tregs_shp %>% filter(Region_Nam == "Mtwara"), fill = NA) +
#   geom_sf(data = mrkta_sf %>% filter(district.region.name == "mtwara"), col = "black", alpha = 0.5, size = 3)+
#   geom_sf(data = mrktCtoAmatch_sf %>% filter(district.region.name.A == "mtwara"), col = "red", shape = 17, size = 1)+
#   # scale_colour_brewer(palette = "Paired")+
#   facet_wrap(~district.region.name)+
#   coord_sf()+
#   theme_bw()
# 
# map_plots[[1]]; map_plots[[2]]; map_plots[[3]]; map_plots[[4]]; map_plots[[5]]; map_plots[[6]]; map_plots[[7]]; 
# map_plots[[8]]; map_plots[[9]]; map_plots[[10]]; map_plots[[11]]

Market Maps (2): Geo-matched Locations

Market A Map 2.1:

sf_use_s2(FALSE)
## Spherical geometry (s2) switched off
mrkta_sf2 <- st_transform(mrkta_sf, crs = st_crs(twards_shp))

twards_shp <- twards_shp %>%
  mutate("District_N" = tolower(gsub(" ", ".", District_N)),
        "Region_Nam" = tolower(gsub(" ",".",Region_Nam)),
        "Ward_Name" = tolower(gsub(" ",".",Ward_Name))
        )

mrktAtwards_match <- mrkta_mapdata %>%
  mutate(twardID = st_nearest_feature(mrkta_sf2,twards_shp)) %>% 
  left_join(
    twards_shp %>%
      mutate(twardID = 1:nrow(.)) %>%
      select(
        twardID,
        "region.name.twards" = "Region_Nam",
        "district.name.twards" = "District_N",
        "ward.name.twards" = "Ward_Name" ,
      ),
    by = c("twardID"))
## although coordinates are longitude/latitude, st_nearest_feature assumes that
## they are planar
# mrktCtwards_match %>% group_by(region.name.twards, district.name.twards) %>% count() %>% View()


mrktAtwards_match_gps <- mrktAtwards_match %>%
  select(
      country,
      region.name.twards,
      district.name.twards,
      ward.name.twards,
      district.region.name.A = district.region.name,
      ward.name.A = ward.name,
      village.A = village,
      market.name.A = market.name.corrected,
      lat,
      long,
      gps.coordinates.of.the.market
    ) %>% 
  mutate(region.district.name.twards = paste(region.name.twards, district.name.twards, sep = "-"))

mrktAtwards_match_sf <- st_as_sf(mrktAtwards_match_gps, coords = c("long","lat"),  crs = 4326)


##################################
# Geo-matched Map, check Regions #
##################################

## Region: id regions with surveyed markets:
regions <- mrktAtwards_match_sf %>% distinct(region.name.twards) %>% pull()
Regions <- str_to_title(regions)

mrktAtwards_match_map <- ggplot() +
  geom_sf(data = tregs_shp, fill = NA, linewidth = 0.1) +

   geom_sf(data = twards_shp %>% filter(Region_Nam %in% regions), aes(fill = Region_Nam), linewidth = 0.1) +
   scale_fill_manual(values=alpha(tanzania_pal,0.7), guide = "none")+
   # geom_sf(data = mrkta_sf, aes(col = district.region.name), alpha = 0.5, size = 3)+
  geom_sf(data = mrktAtwards_match_sf, aes(col = region.name.twards), alpha = 0.5, size = 3)+
  geom_sf(data = mrktAtwards_match_sf, aes(col = district.region.name.A),  shape = 17, size = 1.5)+
   scale_colour_manual(values = rep(market_pal,9))+
   coord_sf()+
   labs(title = "Market Survey A locations: ward-matched (circle) and manual (triangle)",
        fill = "Tanzania region",
        col = "Market region")+
  #theme(legend.position = "none") #+
  theme(legend.position = "bottom", legend.box = "horizontal")
 
mrktAtwards_match_map

#########################################
# Geo-matched Map, check Ward by Region #
#########################################

## Region: id wards with surveyed markets:
wardlist <- mrktAtwards_match %>% select(ward.name.twards) %>% distinct() %>% pull()

map_plots <- list()
for (i in 1:length(regions)) {

  Region_i <- Regions[i]
  region_i <- regions[i]
  
  # For region (i) id wards with surveyed markets according to geo-matched data 
  wardlist_i <- mrktAtwards_match %>% 
    filter(region.name.twards == region_i) %>% 
    select(ward.name.twards) %>% 
    distinct() %>% pull()
  
  map_plots[[i]] <- ggplot() +
    geom_sf(data = twards_shp %>% filter(Region_Nam == region_i), fill = NA) +
    geom_sf(data = twards_shp %>% filter(Region_Nam == region_i,
                                         Ward_Name %in% wardlist_i), 
            aes(fill = Ward_Name, alpha = 0.5)) +

    geom_sf(data = mrktAtwards_match_sf %>% filter(region.name.twards == region_i),
            aes(col = ward.name.twards),  alpha = 0.5, size = 3)+
    # geom_sf(data = mrktAtwards_match_sf %>% filter(region.name.twards == region_i), 
    #         aes(col = ward.name.A),  shape = 17, size = 1.5)+
    # geom_sf(data = mrktCtwards_match_sf%>% filter(region.name.twards == region_i), aes(col = ward.name.3A),  shape = 17, size = 1.5)+
    scale_fill_manual(values = rep(market_pal,10))+    
    scale_colour_manual(values = rep(market_pal,10))+
    facet_wrap(~region.name.twards)+
    coord_sf()+
    theme(legend.position = "bottom")
}

map_plots[[1]]; map_plots[[2]]; map_plots[[3]]; map_plots[[4]]; map_plots[[5]]; map_plots[[6]]; map_plots[[7]];

map_plots[[8]]; map_plots[[9]]; map_plots[[10]]; map_plots[[11]]; map_plots[[12]]; map_plots[[13]]

Market C Map 2.2:

sf_use_s2(FALSE)
mrktc_sf2 <- st_transform(mrktc_sf, crs = st_crs(twards_shp))

twards_shp <- twards_shp %>%
  mutate("District_N" = tolower(gsub(" ", ".", District_N)),
        "Region_Nam" = tolower(gsub(" ",".",Region_Nam)),
        "Ward_Name" = tolower(gsub(" ",".",Ward_Name))
        )

mrktCtwards_match <- mrktc_mapdata %>%
  mutate(twardID = st_nearest_feature(mrktc_sf2,twards_shp)) %>% 
  left_join(
    twards_shp %>%
      mutate(twardID = 1:nrow(.)) %>%
      select(
        twardID,
        "region.name.twards" = "Region_Nam",
        "district.name.twards" = "District_N",
        "ward.name.twards" = "Ward_Name" ,
      ),
    by = c("twardID"))
## although coordinates are longitude/latitude, st_nearest_feature assumes that
## they are planar
# mrktCtwards_match %>% group_by(region.name.twards, district.name.twards) %>% count() %>% View()

mrktCtwards_match_gps <- mrktCtwards_match %>%
  select(
      country,
      region.name.twards,
      district.name.twards,
      ward.name.twards,
      district.region.name.3A = district.region.name,
      ward.name.3A,
      village.3A,
      market.name.3A,
      lat,
      long,
      gps.coordinates.of.the.market
    ) %>% 
  mutate(region.district.name.twards = paste(region.name.twards, district.name.twards, sep = "-"))
  
mrktCtwards_match_sf <- st_as_sf(mrktCtwards_match_gps, coords = c("long","lat"),  crs = 4326)


##################################
# Geo-matched Map, check Regions #
##################################

## Region: id regions with surveyed markets:
regions <- mrktCtwards_match_sf %>% distinct(region.name.twards) %>% pull()
Regions <- str_to_title(regions)

mrktCtwards_match_sf <- st_as_sf(mrktCtwards_match_gps, coords = c("long","lat"),  crs = 4326)

mrktCtwards_match_map <- ggplot() +
  geom_sf(data = tregs_shp, fill = NA, linewidth = 0.1) +

   geom_sf(data = twards_shp %>% filter(Region_Nam %in% regions), aes(fill = Region_Nam), linewidth = 0.1) +
   scale_fill_manual(values=alpha(tanzania_pal,0.7), guide = "none")+
   # geom_sf(data = mrktc_sf, aes(col = district.region.name), alpha = 0.5, size = 3)+
  geom_sf(data = mrktCtwards_match_sf, aes(col = region.name.twards), alpha = 0.5, size = 3)+
  geom_sf(data = mrktCtwards_match_sf, aes(col = district.region.name.3A),  shape = 17, size = 1.5)+
   scale_colour_manual(values = rep(market_pal,9))+
   coord_sf()+
   labs(title = "Market Survey C locations: ward-matched (circle) and manual (triangle)",
        fill = "Tanzania region",
        col = "Market region")+
  #theme(legend.position = "none") #+
  theme(legend.position = "bottom", legend.box = "horizontal")
 
mrktCtwards_match_map

#########################################
# Geo-matched Map, check Ward by Region #
#########################################

## Region: id wards with surveyed markets:
wardlist <- mrktCtwards_match %>% select(ward.name.twards) %>% distinct() %>% pull()

map_plots <- list()
for (i in 1:length(regions)) {

  Region_i <- Regions[i]
  region_i <- regions[i]
  
  # For region (i) id wards with surveyed markets according to geo-matched data 
  wardlist_i <- mrktCtwards_match %>% 
    filter(region.name.twards == region_i) %>% 
    select(ward.name.twards) %>% 
    distinct() %>% pull()
  
  map_plots[[i]] <- ggplot() +
    geom_sf(data = twards_shp %>% filter(Region_Nam == region_i), fill = NA) +
    geom_sf(data = twards_shp %>% filter(Region_Nam == region_i,
                                         Ward_Name %in% wardlist_i), 
            aes(fill = Ward_Name, alpha = 0.5)) +

    geom_sf(data = mrktCtwards_match_sf %>% filter(region.name.twards == region_i), aes(col = ward.name.twards),  alpha = 0.5, size = 3)+
    # geom_sf(data = mrktCtwards_match_sf %>% filter(region.name.twards == region_i), aes(col = ward.name.3A),  shape = 17, size = 1.5)+
  #   geom_sf(data = mrktCtwards_match_sf %>% filter(region.name.twards == region_i)%>% filter(region.name.twards == region_i), aes(col = region.name.twards), alpha = 0.5, size = 3)+
  # geom_sf(data = mrktCtwards_match_sf %>% filter(region.name.twards == region_i), aes(col = district.region.name.3A),  shape = 17, size = 1.5)+
    scale_fill_manual(values = rep(market_pal,10))+    
    scale_colour_manual(values = rep(market_pal,10))+
    facet_wrap(~region.name.twards)+
    coord_sf()+
    theme(legend.position = "bottom")
}

map_plots[[1]]; map_plots[[2]]; map_plots[[3]]; map_plots[[4]]; map_plots[[5]]; map_plots[[6]]; map_plots[[7]];

map_plots[[8]]; map_plots[[9]]; map_plots[[10]]; map_plots[[11]]; map_plots[[12]]; map_plots[[13]]

##################################
# Geo-matched Map, check Regions #
##################################

## Region: id regions with surveyed markets:
regions <- mrktCtwards_match_sf %>% distinct(region.name.twards) %>% pull()
Regions <- str_to_title(regions)

mrktCtwards_match_sf <- st_as_sf(mrktCtwards_match_gps, coords = c("long","lat"),  crs = 4326)

mrktAC_geomatch_map <- ggplot() +
  geom_sf(data = tregs_shp, fill = NA, linewidth = 0.1) +

   geom_sf(data = twards_shp %>% filter(Region_Nam %in% regions), aes(fill = Region_Nam), linewidth = 0.1) +
   scale_fill_manual(values=alpha(tanzania_pal,0.7), guide = "none")+
   # geom_sf(data = mrktc_sf, aes(col = district.region.name), alpha = 0.5, size = 3)+
  geom_sf(data = mrktAtwards_match_sf, aes(col = region.name.twards), alpha = 0.5, size = 3)+
  geom_sf(data = mrktCtwards_match_sf, aes(col = region.name.twards),  shape = 17, size = 1.5)+
   scale_colour_manual(values = rep(market_pal,9))+
   coord_sf()+
   labs(title = "Geo-matched locations: Market A (circle) and Market C (triangle)",
        fill = "Tanzania region",
        col = "Market region")+
  #theme(legend.position = "none") #+
  theme(legend.position = "bottom", legend.box = "horizontal")
 
mrktAC_geomatch_map

##################
# PLOTS BY WARD #
##################

map_plots <- list()
for (i in 1:length(regions)) {

  Region_i <- Regions[i]
  region_i <- regions[i]
  wardlist_i <- mrktCtwards_match %>% filter(region.name.twards == region_i) %>% select(ward.name.twards) %>% distinct() %>% pull()
  
  map_plots[[i]] <- ggplot() +
    geom_sf(data = twards_shp %>% filter(Region_Nam == region_i), fill = NA) +
    geom_sf(data = twards_shp %>% filter(Region_Nam == region_i,
                                         Ward_Name %in% wardlist_i), aes(fill = Ward_Name, alpha = 0.5)) +

    geom_sf(data = mrktAtwards_match_sf%>% filter(region.name.twards == region_i), aes(col = ward.name.twards),  alpha = 0.5, size = 3)+
    geom_sf(data = mrktCtwards_match_sf%>% filter(region.name.twards == region_i), aes(col = ward.name.twards),  shape = 17, size = 1.5)+
    scale_fill_manual(values = rep(market_pal,10))+    
    scale_colour_manual(values = rep(market_pal,10))+
    facet_wrap(~region.name.twards)+
    coord_sf()+
    theme(legend.position = "bottom")
}

map_plots[[1]]; map_plots[[2]]; map_plots[[3]]; map_plots[[4]]; map_plots[[5]]; map_plots[[6]]; map_plots[[7]];

map_plots[[8]]; map_plots[[9]]; map_plots[[10]]; map_plots[[11]]; map_plots[[12]]; map_plots[[13]]

Maps 3: C to geomathced A

mrktCtoAmatch2 <- mrktCtwards_match_gps %>%
  mutate(mrktaID = st_nearest_feature(mrktCtwards_match_sf, mrktAtwards_match_sf)) %>% 
  left_join(
    mrktAtwards_match_gps %>%
      mutate(mrktaID = 1:nrow(.)) %>%
      select(
        mrktaID,
        # district.region.name.A,
        # ward.name.A,
        # village.A,
        # market.name.A,  
        "region.name.twards.A"=region.name.twards,
       "district.name.twards.A"= district.name.twards,
        "ward.name.twards.A"=ward.name.twards
      ),
    by = c("mrktaID"))
## although coordinates are longitude/latitude, st_nearest_feature assumes that
## they are planar
mrktCtoAmatch_gps2 <- mrktCtoAmatch2 %>%
  select(
      country,
      ends_with("twards.A"),
      lat,
      long,
      gps.coordinates.of.the.market
    )
  
mrktCtoAmatch_sf2 <- st_as_sf(mrktCtoAmatch_gps2, coords = c("long","lat"),  crs = 4326)

mrktCAmatched_map2 <- ggplot() +
   geom_sf(data = tregs_shp, aes(fill = Region_Nam), linewidth = 0.1) +
   scale_fill_manual(values=alpha(tanzania_pal,0.7), guide = "none")+
  # market A locations
   geom_sf(data = mrktAtwards_match_sf, aes(col = region.name.twards), alpha = 0.5, size = 3)+
  # market C locations closest to market A locations
   geom_sf(data = mrktCtoAmatch_sf2, aes(col = region.name.twards.A),  shape = 17, size = 1.5)+
    geom_sf_text(data = tregs_shp, aes(label = Region_Nam), size = 2)+

   scale_colour_manual(values = rep(market_pal,9))+
   coord_sf()+
   labs(title = "Market Survey A (circles) and C (triangles) - C matched to geo-A",
        fill = "Tanzania region",
        col = "Market region")+
  theme(legend.position = "right", legend.box = "horizontal")
 
mrktCAmatched_map2
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

ggplot() +
   geom_sf(data = twards_shp %>% filter(Region_Nam %in% c("lindi", "mtwara")), aes(fill = Region_Nam), linewidth = 0.1) +
   scale_fill_manual(values=alpha(tanzania_pal,0.7), guide = "none")+
  # market A locations
   geom_sf(data = mrktAtwards_match_sf%>% filter(region.name.twards %in% c("lindi", "mtwara")), aes(col = region.name.twards), alpha = 0.5, size = 3)+
  # market C locations closest to market A locations
   geom_sf(data = mrktCtoAmatch_sf2%>% filter(region.name.twards.A %in% c("lindi", "mtwara")), aes(col = region.name.twards.A),  shape = 17, size = 1.5)+
   scale_colour_manual(values = rep(market_pal,9))+
   coord_sf(xlim = c(38,39.5), ylim = c(-10,-11.5))+
   labs(title = "Market Survey A (circles) and C (triangles) - C matched to geo-A",
        fill = "Tanzania region",
        col = "Market region")+
  theme(legend.position = "right", legend.box = "horizontal")